A study of the bosonic sector of the two-dimensional 
Hubbard model within a two-pole approximation 
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The charge and spin dynamics of the two-dimensional Hubbard model in the paramagnetic phase 
is first studied by means of the two-pole approximation within the framework of the Composite 
Operator Method. The fully self-consistent scheme requires: no decoupling, the fulfillment of both 
Pauli principle and hydrodynamics constraints, the simultaneous solution of fermionic and bosonic 
sectors and a very rich momentum dependence of the response functions. The temperature and 
momentum dependencies, as well as the dependency on the Coulomb repulsion strength and the 
filling, of the calculated charge and spin susceptibilities and correlation functions are in very good 
agreement with the numerical calculations present in the literature. 

PACS numbers: 71.27, 71.10.f 



I. INTRODUCTION 



The Hubbard model 7 7 7 7 is capable to describe a 
rich variety of behaviors including a wide range of dif- 
ferent spin and charge dynamics 7 . In particular, its 
interactions are thought to be responsible for strong 
antiferromagnetic correlations at half-filling and low 
temperatures 7 . In the presence of doping, the antiferro- 
magnetic correlations remain rather strong although the 
correlation length can get smaller and smaller on increas- 
ing the doping. The possibility of charge order and phase 
separation has also been widely investigated according to 
the fact that one of the mostly used derivative of the Hub- 
bard model, the t-J model, seems to present charge sep- 
aration for a wide range of external parameters 7 . How- 
ever, recent numerical results seems to indicate that the 
two models may have different behavior as far as charge 
correlations are concerned 7 . 

In this manuscript, we first give a fully self-consistent 
treatment of the charge and spin dynamics of the two- 
dimensional Hubbard model in the two-pole approxima- 
tion within the framework of the Composite Operator 
Method (COM) 7 7 . The fermionic and bosonic sectors 
are solved together self-consistently, no decoupling ap- 
proximation is used and the explicit momentum depen- 
dence of the spectra involves third nearest-neighbor sites 
that forces a rather complex and rich momentum depen- 
dence in all physical properties. 

The COM rightfully belongs to the family of the 
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projection methods and is 
based on two main ideas 7 : strongly interacting systems 
should be described in terms of the quasi-particles gen- 
erated by the interactions and the dynamics should be 
bounded to the right Hilbert space through the imposi- 
tion of constraints coming from the Pauli principle. By 
Pauli principle we mean all the relations among opera- 
tors dictated by the algebra 7 . With respect to other 
projection methods the COM has some distinguishable 
peculiarities. In particular, within the COM, there is an 
absolute freedom to choose, as asymptotic fields, those 



that are most suitable with respect to the properties of 
the system we wish to describe. This means that we are 
not bound to any specific recipe to choose them and that 
we can use this freedom to exploit at will the benefits 
coming from one choice or another. We can reproduce 
the results of the other methods in an unique framework 
and also go well beyond. For instance, by choosing suit- 
able asymptotic fields and the closure of their equations 
of motion, we were able to describe the lowest energy 
scale, which is not algebraic in the model parameters, of 
impurity models 7 7 . This result is absolutely precluded 
to other projection methods that uniquely focus on the 
preservation of spectral moments of higher and higher or- 
der. According to this, the method is continuously devel- 
oping as we are constantly seeking, one system after the 
other, both the most suitable asymptotic fields and the 
most effective procedures to determine self-consistently 
the dynamics. 

Once the fermionic propagator is known there are 
several ways to compute the response functions (i.e., 
the retarded propagators of the two-particle excitations: 
charge, spin, pair, . . . ). Many techniques are related to 
the possible diagrammatic expansions of the two-particle 
propagators in terms of the single-particle one (i.e., the 
fermionic propagator). However, when operators with 
non-canonical commutations are involved the only fea- 
sible approach is based on the one-loop approximation. 
The complicated algebra of the composite operators in- 
validates the Wick theorem and, consequently, docs not 
allow any simple extension of decoupling schemes and 
more involved diagrammatic approximations 7 7 . Ac- 
cording to this, we have developed and widely applied 
a standard procedure to use, by means of the equations 
of motion approach, the one-loop approximation for com- 
posite operators 7 . 

In this manuscript, we consider another way to tackle 
the problem: the two-particle excitations can be consid- 
ered as a new sector in the dynamics of the system and we 
can choose also for them a suitable asymptotic basis alike 
it has been done for the fermions. This gives a new set of 
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equations obeyed by the two-particle Green's functions 
and the appearance of zero-frequency constants and un- 
known correlators. Also in this case, the enforcement of 
the constraints deriving from the Pauli principle allows to 
compute all the parameters and to fix the representation 
of the Hilbert space 7 . 

Within the framework of the COM, both methods have 
advantages and disadvantages. The one-loop approxima- 
tion becomes exact in the non-interacting limit, well de- 
scribes the incoherent behavior of the two-particle prop- 
agators and establishes a tight connection between the 
one- and two- particle propagators. These are really rel- 
evant properties once we wish to describe the bosonic 
excitations starting from its fundamental constituents: 
the electrons. The Fermi surface bending and nesting 
and the position of the van Hove singularities strongly 
influence the charge and spin response functions. Ac- 
cording to this, we managed to give an explanation for 
the spin magnetic incommensurability issue 7 and the 
overdoped transition of the cuprate superconductors 7 . 
On the other hand, the one-loop approximation is not 
adequate to describe the system in the proximity of or- 
dered phases as the dynamics of the bosonic excitations 
is mainly described in terms of scattering of elementary 
electronic excitations. This practically induces so strong 
finite life-time effects to prevent any possible softening 
of the bosonic modes. As discussed in the above, any 
possible extension involves so complicated diagrammatic 
expansions to be practically unfeasible. 

As regards the pole approximation for the two-particle 
propagators we have obvious advantages like: the pos- 
sibility to easily get the spectra and the analytical ex- 
pressions of correlation functions and susceptibilities; the 
capability to study instabilities (i.e., the softening of the 
modes) in the whole range of model and external pa- 
rameters; the possibility to consider the bosonic excita- 
tions as the media of new interactions among the elec- 
trons. In this paper, we show that it is possible to get 
spin antifcrromagnetic correlations and weak charge or- 
dering tendency at commensurate filling in exceptionally 
good agreement with the numerical results present in the 
literature. On the other hand, the pole approximation 
is based on a description of the bosonic excitations as 
true quasi-particles: the two-particle properties are en- 
tirely controlled by the dynamics, which is only weakly 
renormalized by the fermions; the single-particle proper- 
ties (e.g., Fermi surface shape, position of the van Hove 
singularity, . . . ) do not influence significatively the re- 
sponse function behaviors; the finite life-time effects are 
completely neglected and the tendency towards ordering 
(i.e., softening) is sometime exaggerated. Anyway, the 
use of the Green's function formalism for the bosonic 
sector opens the possibility to explore another really rel- 
evant issue: the ergodicity of the bosonic dynamics and 
the presence of zero-frequency constants in the expres- 
sion of the casual Green's function and of the correlation 
functions 7 7 7 . In this manuscript, we decided not to 
pursue this analysis and to fix the zero-frequency con- 



stant values by means of ergodicity requirements in ac- 
cordance with the general understanding of bulk systems. 

As we can see, the two methods are effectively comple- 
mentary and can be used to describe the spin and charge 
dynamics of the system in different region of the param- 
eter space according to the relevance and prevalence of 
localization and ordering (two-pole) with respect to de- 
localization and symmetry (one- loop). 

It is also worth noting that the pole approxima- 
tion allows, at least in principle, to get a completely 
self-consistent solution putting together fermionic, spin, 
charge and pair propagators 7 . The Pauli principle could 
be then used to get also the zero-frequency constants in 
self-consistency and definitely fix the Hilbert space, as 
described in Ref. ? . 



II. THE HUBBARD MODEL AND THE 
FERMIONIC SECTOR 

The Hubbard model is described by the following 
Hamiltonian 

h = E (*u - n *y) ct « c V) + u E "t (0 n i (<) 

ij i 

(2.1) 

where (i) = ^cj (i) cj is the electronic creation 

operator in spinorial notation at the site i [i = (i, t)] and 
n a (i) = c\ it) c a (i) is the number operator for spin a at 
the site i; fi is the chemical potential and U is the on-site 
Coulomb repulsion. 

The matrix ty describes the nearest-neighbor hopping; 
in the 2D case we have ty = — 4tay, where 



a ij = ^£ eik(i " j M k ) 



(2.2) 



is the projector on the nearest-neighbor sites and a (k) = 
i [cos (k x a) + cos (k y a)) and a is the lattice parameter. 



We choose the following fermionic basis - ■ 1 



(2.3) 



where £ (i) = [1 — n (i)] c (i) and r){i) = n (i) c (i) are the 
Hubbard operators. (i) satisfies the following equation 
of motion 



-Atc a (i) — 4* 7r (i) 
— (A* — I7)»7 (*) + 4t7T (*) 



(2.4) 



where c 7 (i, t) = ^\ 7y c (j, t) [7y stands for any projector 
on the j neighbor sites of i; see Appendix] and tt (i) = 
\o» TV (i) c a (»)+c (i) [ct« (») c (t)] . n„(i) = c t(i) <t„ c(i) 
are the charge (/x = 0) and spin (/i = 1,2,3 ) density 
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operators, with ay, = (1, ct), ct^ 1 = (—1,0*) and a are the 
Pauli matrices. 

Let us project the source J (i) on the chosen basis 



J(i,t)e*5^e(i,j) *(j,t) 



(2.5) 



Accordingly, the energy matrix e (i, j) is defined through 
the equation 



rn 



(i,j) = ^ £ (i,l)7(l,j) 



(2.6) 



where the m-matrix and the normalization matrix I have 
the following definitions 



m(i,j) = ({j(i,t),*t(j,t)}) 

j(ij) = ({*(i,t),*ta*)» 



(2.7) 
(2.8) 



It is worth pointing out that in Eq. Ij2.7|l J(i) is the to- 
tal current given in Eq. I|2.4|l and not the approximated 
one. After Eq. (|2.5() . the Fourier transform of the ther- 
mal single-particle retarded Green's function G(i,j) — 
(R [ty (i) ^ (j)] } satisfies the following equation 



[w-e(k)]G(k,«) = I(k) 



(2.9) 



9 V f ? 



The straightforward application of this scheme' 
gives that, in the paramagnetic phase, I (k) has diagonal 
form with J n = 1 — n/2 and I 2 2 = n/2 ((n a (i)) = ^) 
and that the m-matrix depends on three parameters: the 
chemical potential fj, and two correlators 

A = (C (i) e (i)> - (V a (*) ^ (0) (2-10) 

P = W (»')> - (t c T (0 c i ( J )] Q 4 (0 4 (*)> 

(2.11) 

The three parameters /x, A and p are functions of the 
external parameters n, T and Z7 and can be determined 
self-consistently through a set of three coupled equations 



n = 2 [1 - (c (i) 



«(<) W> 



2+ 
(<))• 





(0>] 

- (?7 Q (*) *? f (<)> 



(2.12) 



The first equation fixes the particle number, the second 
one comes from the definition of A and the third one as- 
sures that the solution respects the Pauli principle (i.e., 
the local algebra)' . In this latter equation resides the 
main difference with all the other two-pole approxima- 
tions. This equation: allows to fix the representation 7 ; 
implements the particle-hole symmetry in the solution 7 ; 
avoids uncontrolled decoupling or further approximations 
on higher order correlators. Using this procedure is pos- 
sible to find two solutions: one with a p positive and of 
order of the filling n and another with p manly negative 
and rather small. The main difference between these two 



solutions resides in the strength of the antiferromagnetic 
correlations 7 7 . 

It is worth noting that this set of coupled self- 
consistent equations gives the exact solution in the 
atomic and in the non-interacting cases as well as for 
the interacting two-site system 7 . According to this, we 
are confident to reproduce at least the two most relevant 
scale of energies in the system: the Coulomb repulsion 
U and the exchange energy J. The latter is already well 
defined on the two-site system that is the minimal cluster 
where J becomes effective. 

Within this calculation scheme, the fermionic solution 
is known in a fully self-consistent manner and without 
opening the bosonic sector. Once we have the electronic 
Green's function we can get all single-particle, local and 
thermodynamic properties straightforwardly. In the last 
years, by means of the COM in the above described ap- 
proximation, we got results in excellent agreement with 
numerical and exact solutions as regards many lattice 
and impurity systems 7 



III. CHARGE AND SPIN RESPONSE 
PROPERTIES 

As stated in the introduction we choose to compute the 
charge and spin response functions by studying the causal 
Green's function 7 G^ — (T [n^ (i) n M (j)]). As we 
widely discussed in Ref. ? , to obtain a correct descrip- 
tion of the bosonic properties is necessary to compute 
first the causal Green's function and then derive from this 
latter all other propagators and correlators. The reason 
of this lies in the fact that the zero-frequency constants 
do not explicitly contribute to the retarded functions, 
although there is an implicit dependence through the 
self-consistent determination of the internal parameters. 
Starting from the retarded function would lead to ignore 
the zero-frequency constants and will give wrong re- 
sults. Once we know the Fourier transform of G^*) 
that is G^> (k, u), we can find spin and charge suscep- 
tibilities x (k, w) = —F (R [rip (i) n M (j)}} and correla- 
tion functions C^ (k, u>) = F (n^ (i) (j)) by means of 
the well-known formulas 



>t X (M) (k,w) 
X M (k,w)' 
C M (k,w) = 



= -8t \G M (k,uj) 



= — tanh — 3 
2T 



1 



- tanh — ) 
27V 



G^ (k,w 
G M (k,w) 



(3.1) 
(3.2) 
(3.3) 



where T, R and F are the causal and retarded operators 
and the Fourier transform, respectively. 

As widely discussed in the introduction, in this 
manuscript we will study the spin and charge channels 
of the bosonic sector by using a pole approximation. Let 
us write the equation of motion for the operator (i) 



. d 



(*) = -U Pli (<) 



(3.4) 
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where 

p» W = ct (*) ^ c " W - cW W °M c W 



(3.5) 



The bosonic basis has to be chosen in order to be com- 
patible with the fcrmionic one and with a non-local com- 
ponent as we wish to take into account, at least par- 
tially, the derealization driven by the kinetic term of the 
Hamiltonian. According to this, the most natural choice 
is a two-component basis and, in particular, that directly 
generated by the hierarchy of the equations of motion. 
This will assure that the first four bosonic spectral mo- 
ments have the correct functional form 7 . Therefore, we 
take as bosonic basis the following one 



(*) 



(3.6) 



The equation of motion of (i) is the following one 
. d 



Ulp(i) + UKp(i) (3.7) 

where the higher-order bosonic operators are defined by 

k m (i) = c f (i) cr M ij a (i) - (i) cr M c a (i) 

+^ a {i)a^c(i)-c^(i)a^(i) (3.8) 

J„(i) = c+(i) <r M c" 2 (i) + C t« 2 (i) a„c(i) 

-2ct Q (i) ^ c a (») (3.9) 

and the definition of c" 2 (i) can be found in Appendix. 

Using the same procedure used for the fermions, we 
have 



d 



(3.10) 



rri 



(m) 



(k) 



roft> (k) 







22^ (k) 



where 

I<£>(k) = 4[l-a(k)]C° 



(k) 

» 



-4t/^ (k) 



m^(k) = -4tJ, M „„ (k) + (k) 



(3.15) 

(3.16) 
(3.17) 
(3.18) 



The exact expressions of the normalization matrix entries 
and the definition of the parameters they depend on can 
be found in the Appendix. The energy matrix (k) 
has off-diagonal form with non-zero elements 



etf (k) = -4t 



-21 



(k) 



» 



(k) 



I® (k) 



(3.19) 
(3.20) 



For the sake of simplicity, we will now extend the pre- 
vious used notation for the bosonic causal Green's func- 
tion = (T [n,,. (i) n M (j)]) to the complete 2x2 
matricial one, that is, is hereafter defined as 
(T [Nfj, (i) (j)])- We will also use the accordingly ex- 
tended notation for the correlation function 
They have then the following expressions 



G<"> (k.w) = (k) 8(lo)T^ 



+ ° M (k) ■ 



1 + /bH 



71=1 



w-w&° (k)+i<5 



-jV".") (k) 



/bH 



n=l 



(k) -iS 



(3.21) 



where (i, j) is given by 

ro<"> (i,j) = ^ £ W (i,l) (l,j) 



(3.11) 



and the normalization matrix and the m^-matrix 
have the following definitions 



(/') 



(U) 



i^JMM),jvt(j,f) 



(3.12) 
(3.13) 



As it can be easily verified, in the paramagnetic phase the 
normalization matrix does not depend on the index 
/U; charge and spin operators have the same weight. The 
two matrices 1^ and have the following form in 

momentum space 7 



(k) 




(3.14) 



C M ( k)W ) = ^!5( 2 ) (k) ^r, 



2 



+ 2tt^ [w - w<"> (k)] [1 + M">)] ° M (k) 

71=1 



(3.22) 

where r M is the zero- frequency constant 7 and /b(w) = 
-;rs — r is the Bose-Einstein distribution function. In this 



manuscript, we will use the ergodic value (i.e., r llM = 
Sfj,o n 2 ) for the zero- frequency constant as explained in 
the introduction. The energy spectra are given by 

(3.23) 
(3.24) 



(k) = (-)"<>> (k) 



(/<■) 



(k) = V^ 



12 (k) £21 



(m) 



(k) 



and the spectral functions have the following expression 



>(k) 



(k) = ^ 



(k) 



; (k) 



4? po 

3? (k) 



(3.25) 
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As it can be seen from the expressions given in Appendix, 
the Green's function and the correlation function depend 
on various parameters, static correlation functions, that 
must be self-consistently calculated. A subset of param- 
eters, C Q , C A , C^, E 13 and E v , are of fermionic nature 
and can be computed through the fermionic Green's func- 
tion. The negative p solution will be used in order to get 
enhanced antiferromagnetic correlations. The remaining 
parameters, a^, b^, and d^, are of bosonic nature, 
but they cannot be expressed in terms of the bosonic 
Green's function under study as they belong to higher 
order propagators. As in the fermionic sector, we can 
avoid studying complicated higher order propagators re- 
quiring the fulfillment of the Pauli principle and of other 
symmetry requirements. Four equations will be used to 
fix these parameters: one equation comes from the Pauli 
principle and other three from the general properties of 
the bosonic spectra a>„ (k) for small momenta (i.e., for 
k — > where k = « /k% + ky). The Pauli principle 7 gives 



(n(i) n(i)) 
(n k (i) n k (i)) 



2D 
2D 



k = 1, 2, 3 



(3.26a) 
(3.26b) 



where D = (n^ (i) nj, (i)) = § — (n (i) (i)) is the double 
occupancy. From the continuity equation 7 it follows that 
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lini iu#> (k) = #> k s (3.27) 

k — >U 

where s > 1 and cifi is the velocity. Let us analyze the 
expression for u>„ (k) . The function ro^' (k) at small k 
can be cast in the following form 

m$ (k) = m£° + {kaf + (k a) 4 

+m^ ) (k a) 4 sin 2 (2(f> k ) + 0((k a) 6 ) (3.28) 

where <pk = arctan ^ . The function F$ (k) behaves as 
(ka) 2 C a at small k. Therefore, to satisfy the continuity 
equation we must put 

= = (3.29) 

Moreover, as the susceptibility has to be single-value 
at k = we have also to require that m| = 0. The 
coefficients of (k) in the limit of small k have the 
following expressions (see Appendix) 

= U(-a^ + i& M + i Cp + 1^) (3.30) 

m { f = j (2a M -c li -d li -2D- 2E 7 ') (3.31) 
mf = ~\t{C a -2C" 1 + C X ) 

+ ^-(a„ + c, Jl -2d fl -D + 6EP-7E*) (3.32) 



FIG. 1: a3 and ao as functions of the filling n for T = and 
U = 4 and 8. 

According to this, we can write the following algebraic 
relations for the parameters 6^, c M and d M 

b,, = a M + 3D + E v + 2E (i 

~ 6 JJ {° a + C A - 2C"*) (3.33) 
c M = a^-D + E' 1 - 2E 13 

+ 6 JJ i " + C A - 2C" 1 ) (3.34) 
d^ = a^-D- 3E r > + 2E 

(C a + C x - 2C") (3.35) 

and use the Eq. (|3.26|l to compute the parameter 
self-consistently. In Fig. ^ we report the behavior of 03 
and ao as functions of the filling n for T = and U = 4 
and 8. The behavior of 03 reveals a strong dependence 
on both filling and Coulomb repulsion of the intensity of 
spin correlations. In particular, at half-filling we have the 
maximum dependence on U. ao, instead, is practically 
featureless except for a region near half filling, whose 
extension is controlled by the strength of the Coulomb 
repulsion, where rapidly and enormously increases with 
a slope that again depends on U. This latter behavior 
results in a strong enhancement of the charge correlations 
in the proximity of the Mott-Hubbard metal-insulator 
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FIG. 2: The spin o/ 3 '(k) and charge ^"'(k) spectra as func- 
tions of the momentum (k x = k y ) for n = 1, 0.9, 0.8, U — 4 
and 8 and T = 0; the qMC + data (10 x 10) for w (3) (k) at 
U = 4, n = land T — are from ? . 



FIG. 3: The uniform static spin susceptibility xo as function 
of the temperature T for U = 4, n = 1, 0.75 and 0.25; the 
qMC data (8 x 8) are from ? . 




(0,0) («.*) k M) (0,0) 



transition. 

It is necessary to report that this analysis can be con- 
sidered an extension and a completion of that done in 
Ref. ? . The main differences are related to the use of 
causal propagator in place of the retarded one and to the 
exploitation of the hydrodynamics constraints to fix the 
parameters coming in the energy spectra whenever we 
wish to retain the complete dependence on the momen- 
tum. 



IV. RESULTS 
A. Spin and charge spectra 

The spin and charge spectra, as functions of the mo- 
mentum, are reported in Fig- Elfor n = 1, 0.9, 0.8, U = 4 
and 8 and T = 0. As regards the spin spectrum, COM 
result is in fair agreement with the quantum Monte Carlo 
data 7 (10 x 10) except for k = (tt, n) = Q. The very 
small value reported by the numerical data at Q can be 
understood as a consequence of overestimated antiferro- 
magnetic correlations (i.e., the antiferromagnetic correla- 
tion length actually exceeds the cluster size, see Fig. [jJJ. 
COM results, instead, are obtained with paramagnetic 



FIG. 4: The spin susceptibility x(k) as function of the mo- 
mentum for U = 4, T = and n = 1, 0.9 and 0.8. 



boundary conditions. The minimum at Q in the spin 
spectrum is the clearest possible evidence that we have 
quite strong antiferromagnetic correlations in our solu- 
tion. The doping is quite efficient in reducing the in- 
tensity of them. On the contrary, they significatively 
increase with the Coulomb repulsion according to the 
stronger and stronger influence of the exchange energy J 
in the strongly interacting regime. The charge spectrum 
shows an enhancement of the velocity with decreasing 
doping and increasing Coulomb repulsion, that is, in the 
proximity of a Mott-Hubbard metal-insulator transition, 
which would have as signature the divergency of the for- 
mer. 



B. Spin susceptibility 

The dynamical spin susceptibility Xs (k, to) can be ob- 
tained by Eqs. (|3.1fl and l|3.2|) with /i = 3 and has the 
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FIG. 5: The spin correlation function S(k) as function of the 
momentum for (7 = 4, T = 0.2 and n = 0.8, 0.33 and 0.19; 
the qMC data (8 x 8) are from ? . 



FIG. 6: The spin correlation function S(k) as function of the 
momentum for U = 8, T = 0.2 and T = 0.57 and n = 1, 0.45 
and 0.2; the qMC data (8 x 8) are from ? . 



expression 



Xs (k,w) = - 



2 



r (™,3) 
'11 



(k) 



w - wi 3) (k) + i 5 



(4.1) 



According to this, the static x (k) = lim w ^ Xs (k, oj) and 
the static and uniform xo — linik^o X (k) spin suscepti- 
bility are given by 



X(k) 



Xo 



{4[l-a(k)]C°r 



(4C a ) 



a^2 



24t(C« - C^ 1 ) - [7(c 3 + 4EP) 



(4.2) 
(4.3) 



X0j a s a function of the temperature, is reported in Fig. [2 
for U = 4 and n = 0.25, 0.75 and 1. COM results 
are in very good agreement with the quantum Monte 
Carlo ones 7 (8 x 8) for n = 0.25 and 0.75. For n = 1 
and low temperatures, our paramagnetic solution can- 
not reproduce the overestimated antiferromagnetic cor- 
relations present in the numerical results. Anyway, our 
spin susceptibility \ (k) and our spin correlation function 
S (k) (see next section) present a large enhancement at 
Q on reducing the doping (see Fig.0} and increasing the 
Coulomb repulsion (see Figs. and 0) showing that also 
COM results present well developed antiferromagnetic 
correlations although they should be compatible with the 
chosen paramagnetic solution. It is worth noting that the 
presented results are in better agreement with quantum 
Monte Carlo data than the random phase approximation 
and the COM within the one-loop approximation (see 
Ref. ? and references therein). 



C. Spin correlation function 

The spin correlation function is defined as 

S (i, j) = (n 3 (i, t) n 3 (j, i)> = F- 1 [S (k)] (4.4) 
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FIG. 7: The spin correlation function S(k) as function of the 
momentum for n — 0.875, (7 = 8 and T = and 0.17; the 
Lanczos data (4x4) at T = are from ? . 



where F 1 stands for the Fourier anti-transform and the 
structure factor reads as 



S(k) = - 



2tlg> (k) 
wW (k) 



coth ■ 



,(3) 



(k) 



2T 



(4.5) 



The behavior of S (k) , as function of the momentum, is 
reported in Figs. [SI El arid El in comparison with some nu- 
merical data 7 7 7 for different values of filling, Coulomb 
repulsion and temperature. We have a very good agree- 
ment with the numerical results for sufficiently high val- 
ues of the doping for all shown values of the Coulomb 
repulsion. In the proximity of half-filling the numerical 
data suffer from a saturation of the antiferromagnetic 
correlation length 7 that becomes comparable with the 
size of the cluster. For U — 4 and n — 0.8 (see Fig. |5J|, 
the correlation length is slightly smaller than the size 
of the cluster: our solution results capable to describe 
this situation fairly well (the height of the peak at Q 
is exactly reproduced and the momentum dependence is 
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FIG. 8: The spin correlation function at Q (S'(Q)) as function 
of the inverse of U for n = 0.875 and T — and 0.17; the 
Lanczos data (4 x 4) at T = are from ? . 



qualitatively correct, again practically exact along the di- 
agonal) except for the exact value of the numerical data 
along the main axes. This is probably due again to an 
overestimation of the correlations by the numerical anal- 
ysis owing to the finite size of the cluster. For U = 8 
and n = 1 (see Fig. [5} and 0.875 (see Fig. 0, in order 
to reproduce the numerical data we need to increase the 
temperature as to decrease our value of the correlation 
length and match that of the numerical analysis, which is 
stuck at the saturation value due to the finiteness of the 
clusters. The results of such a procedure are astonishing, 
we manage to exactly reproduce the numerical data for 
any value of the momentum, and not only at Q, reveal- 
ing the correctness and power of our approach and the 
limitations of the numerical analysis. According to this, 
we wish to point out that the numerical data need to be 
carefully and cleverly interpreted in order to obtain from 
them sensible and effective information. The behavior of 
the spin correlation function as a function of 1/U oc J 
(the exchange energy) is shown in Fig. EI Again, in order 
to obtain results comparable with the numerical ones 7 
we need to use an higher temperature: at T — and for 
high enough values of U, our results show a divergency in 
the correlation length that is extraneous to the physics of 
a finite system. By using the same temperature of Fig. [7| 
(the Lanczos data have the same source), we get a very 
good agreement for any regime of interaction: our solu- 
tion properly describes also the low-energy dynamics of 
the spin system. 

In Fig. we report the behavior of S(i,j) along the 
three principal directions in the lattice for U = 4, T = 0.1 
and (top) n = 1 [(bottom) n = 0.5]. The quantum 
Monte Carlo results 7 at n = 1 present an antiferromag- 
netic correlation length as large as the size of the cluster. 
The correlation along the principal axes [(0,0) — * (i x ,Q) 
and (5,0) — > (5,i y )] is antiferromagnetic and is ferro- 
magnetic along the diagonals [(0, 0) — > (i, i)] as in an or- 
dinary two-dimensional Need phase. COM results show 
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FIG. 9: The spin correlation function S(i,j) along the prin- 
cipal directions for U = 4, T = 0.1 and (top) n = 1 [(bottom) 
n = 0.5]; the qMC data (10 x 10) are from ? . 



exactly the same behavior although the correlation length 
is much smaller: we analyze the paramagnetic phase and 
for U = 4 we still not have so well developed antifer- 
romagnetic correlations. The on-set of an antiferromag- 
netic phase (i.e., to have an antiferromagnetic correlation 
length as large as the size of the system) for a finite sys- 
tem seems possible for any finite value of U at half-filling, 
while, for an infinite system, it could be related to the 
existence of a critical value of the interaction U that, in 
our case, falls between 4 and 8. Actually, our study of the 
antiferromagnetic phase 7 confirm that our critical value 
is C/ c = 6. At n = 0.5 the agreement becomes quantita- 
tive as the strong antiferromagnetic correlations present 
at half filling completely disappear. 



D. Charge correlation function 

The charge correlation function is defined as 

A^(i,j) = (n(i,t)n(j,t))=^ 1 [^ (k)] 
where N (k) reads as 

(k) 



2T 



(4.6) 



(4.7) 
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FIG. 10: The charge correlation lunction iV(k) as a function 
of the momentum for U = 8, T = 0.125 and 0.25 and n — 
0.155, 0.2 and 0.5; the qMC data (8 x 8, 12 x 12, 16 x 16) are 
from ? . 
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FIG. 11: The charge correlation function N(r) as a function 
of the distance for n = 0.2, U = 8 and T = 0.25; the qMC 
data (16 x 16) are from ? . 



N (k) is reported in Fig.^J as a function of the momen- 
tum, for various fillings and temperatures and U = 8. We 
have again a very good agreement with quantum Monte 
Carlo results 7 for all shown values of the external pa- 
rameters and of the momentum. The enhancement at 
k = Q = M for n = 0.5 can be interpreted as the man- 
ifestation of a rather weak ordering of the charge with 
a checkerboard pattern. COM results manage to repro- 
duce such double peak structure showing a capability to 
quantitatively describe, also in a translational invariant 
phase, rather strong charge correlations. 

In Figs.lllland ll2l we report the behavior of N (r), as 
a function of the distance r = \/ i 2 + j 2 , for U — 4 and 
12, T = 0.01 and n = 8/9 and U = 8, T = 0.25 and n = 
0.2, respectively. COM results are in good quantitative 
agreement with the numerical results' " showing once 
more that the charge dynamics is really well described by 
our solution. In Fig. 1131 N (i, i ax ) is shown as a function 
of the Coulomb repulsion U for n = 8/9 and T = 0. The 
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FIG. 12: The charge correlation function N(r) as a function 
of the distance for n = 8/9, [7 = 4 (top) [U = 12 (bottom)] 
and T = 0; the Lanczos data (4 x 4) are from ? . 




FIG. 13: The charge correlation function N (i,i ax ) as a func- 
tion of U for n = 8/9 and T = 0; the Lanczos data (4 x 4) 
are from ? . 



agreement with Lanczos data 7 is quite good and gets 
better and better as U increases. 



V. CONCLUSIONS 

An analytical description of the charge and spin dy- 
namics of the two-dimensional Hubbard model in the 
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paramagnetic phase has been presented within a two- 
pole approximation in the framework of the COM. The 
hydrodynamics constraints as well as the Pauli princi- 
ple requirements have been embedded in the fully self- 
consistent solution by the very beginning and any de- 
coupling has been avoided. The antiferromagnetic cor- 
relations are really well described together with some 
weak charge ordering tendency at commensurate filling. 
Spin spectrum, static uniform spin susceptibility, spin 
and charge correlation functions are in very good agree- 
ment with the numerical results present in the literature 
and clearly state the reliability of the proposed proce- 
dure. 
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a M = 2< C t(i) a^c a {i) c + (») ^c a (i)) 

-(c Qt (i) a M a x a^c a (i) n x (i)) 
b, t = 2(ct(i) CT/iC t(i) <7„[c(i) c(i)] a ) 
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where we used the notation 

i Q = (i x + a,i y ,£) 
i'' = (i x + 2a, ij,, i) 
i' 3 = (i x + a,i y + a, t) 
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APPENDIX 

We have the following expressions for the m^-matrix 
entries 



Kp. ( k ) 



3 r 
4^ 



: (k)] (12C Q + C x + 6(7") 



(k) 



--[l-ry(k)](C Q + C A + 2C") 

+ i[l-A(k)]C A + ^[l- M (k)]C" 

-3[l-/3(k)] (C Q + C") 

-2 [1 - a (k)] D 

+ [1 - 2a (k)] (2£^ + £") 

(k) E 71 + 2/3 (k) £^ 
+ [1 - 2a (k)] a p 

+i[6 M + 2/3(k) c^ + r7(k) d M ] 



The following definitions have been used 

C a = (c a (i)cHi)) 

C x = (c A (7)^(7)) 

C" = <c" (i) c f (»)) 

= (c^ (z) 77+ (0) 

= (c 17 (i) 77 f (»)) 



(1) 



(2) 

(3) 
(4) 
(5) 
(6) 
(7) 



The functions /3y, 77y, fj,ij and A^, the projectors on 
the second, third, fourth and fifth nearest neighbors, re- 
spectively, have the following expressions in momentum 
space 



/?(k) = 


- {cos [a 4 


- ky)] + cos [ 


a (fcr - 


fcv)]} (16) 


77 (k) = 


- [003(20/03;) - 


f cos (2a k y ] 




(17) 


M(k) = 


— {cos [a (2^ 


+ k y )]+ COS 


[a (k x - 


f 2^)] 




+ cos [a (2fe x - 


- ky)] + cos 


a (fej. - 


-2fc„)]} (18) 


A(k) = 


- [cos^afe^) • 


f cos (3a &j,) 


'] 


(19) 



The following relations hold 



C° 2 « = i[ C (7)+2^(7)+ C "(z)] 

c" 3 (i) = ^ [9c Q (7) + c A (7) + 6c" (7)] 



(20) 
(21) 
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